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ABSTRACT 


Genomic analyses of bacterial isolates are necessary to monitor the prevalence of antibiotic 
resistance genes and virulence determinants. Herein, we provide a comprehensive genomic 
description of a collection of 339 Staphylococcus aureus strains isolated from patients with 
bacteremia between 2014 and 2022. Nosocomial acquisition accounted for 56.6% of episodes, 
with vascular catheters being the predominant source of infection (31.8%). Cases of fatality 
(27.4%), persistent bacteremia (19.5%) and diagnosis of septic emboli (24.2%) were 
documented. During the COVID-19 pandemic, we observed a 140% increase of the episodes of 
S. aureus bacteremia per year, with a concomitant increase of the cases from nosocomial origin. 
This prompted us to investigate the existence of genetic features associated with S. aureus 
isolates from the COVID-19 pandemic. While genes conferring resistance to B-lactams (blal-blaR- 
blaZ), macrolides (ermA, ermC, ermT, mphC, msrA) and aminoglycosides (ant(4’)-la, ant(9)-la, 
aph(3’)-Illa, aph(2")-Ih) were prevalent in our collection, detection of the msrA and mphC genes 
increased significantly in pandemic S. aureus isolates. Similarly, we observed a higher prevalence 
of isolates carrying the genes encoding the Clumping Factors A and B, involved in fibrinogen 
binding. Of note, macrolides were extensively used as accessory therapy for COVID-19 and 
fibrinogen levels were usually elevated upon SARS-CoV-2 infection. Therefore, our results reveal 
a remarkable adaptation of the S. aureus isolates to the COVID-19 pandemic context and 
demonstrates the potential of whole-genome sequencing to conduct molecular epidemiology 


studies. 
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INTRODUCTION 


Bacteremia is a life-threatening condition with high morbidity and mortality rates that can elicit 
a systemic host response known as sepsis (Viscoli 2016). These episodes often result in 
complicated bacteremia, usually defined by the presence of attributable mortality, the 
development of hematogenous embolisms (i.e. endocarditis, discitis, osteomyelitis of systemic 
abscesses) or the persistence of viable bacteria in blood after three or more days after proper 
antibiotic treatment (Gudiol et al. 2015). Staphylococcus aureus is among the top species causing 
bacteremia-associated mortality (Kern and Rieg 2020) and persistent bacteremia (Wiggers et al. 
2016), which in turn is linked to a higher risk of metastatic spread (Khatib et al. 2006). There is 
increasing awareness that patients who survive sepsis often have long-term physical, 
psychological and cognitive disabilities with significant health care and social implications 


(Singer et al. 2016). 


Several host factors such as certain comorbidities, the setting of the infection onset and the 
severity of sepsis, have been associated with the high mortality rates observed in S. aureus 
bacteremia (van Hal et al. 2012; Yilmaz et al. 2016). However, the patient’s age is the most 
consistent predictor of mortality, with older patients being twice as likely to die (van Hal et al. 
2012). Also, certain strain characteristics such as methicillin resistance have also been shown to 
impact the clinical outcome (Lee et al. 2018). In contrast, the role of specific genetic traits of S. 
aureus in the evolution of bacteremia is much less understood. Recent studies suggest that S. 
aureus genetic factors predictive of infection evolution may be lineage specific, which would 


hinder their identification (Gasch et al. 2013; Recker et al. 2017). 


S. aureus has also been described as a leading cause of secondary infection during past viral 
pandemics, significantly increasing patient mortality rates (Chertow and Memoli 2013). In this 
regard, recent reports have noted a higher incidence of bacteremia by S. aureus at the beginning 
of the COVID-19 pandemic (Falces-Romero et al. 2023). In the same line, our group has recently 
reported a rise of vascular catheter-related bacteremia in early 2020. Moreover, we also 
observed a significant increase in cases of S. aureus bacteremia (Gasch et al. 2022). However, 
while several studies have addressed the incidence, prevalence and clinical outcomes of SARS- 
CoV-2 and S. aureus co-infection (Adler et al. 2020; Adalbert et al. 2021), the genomic 
differences between such strains and those isolated in the pre-pandemic times have not been 


explored. 


Here, we investigated the clinical aspects, epidemiology and genomic features of 339 S. aureus 


strains isolated from patients with bacteremia at the Parc Taulí University Hospital before and 
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87 during the COVID-19 pandemic (2014-2022). This analysis allowed us to shed light into significant 
88 associations between S. aureus lineages, antibiotic resistance genes, virulence factors and the 
89 clinical outcome of bacteremia. Comparative genomics analysis pointed out that some virulence 


90 factors and antibiotic resistance genes were enriched in pandemic S. aureus isolates. 
91 

92 MATERIAL AND METHODS 

93 Cohort and strain isolation 


94 Patients and setting. All consecutive adult patients (218 years) diagnosed with S. aureus 

95 bloodstream infection from July 2014 to December 2022 at the Parc Tauli University Hospital 

96 (Sabadell, Spain) were retrospectively included in this study. 

97 Clinical data. Clinical information was recorded including patient demographics, comorbidities 

98 and other clinical characteristics, source of bacteremia, place of acquisition, severity of sepsis, 

99 antibiotic therapy and infection outcomes. The modified Charlson score (Charlson et al. 1987) 
100 was used to assess patient morbidities. Setting of acquisition was classified according to 
101 modified Friedman criteria (Friedman et al. 2002) as: community, in-hospital acquired or 
102 healthcare-related. 30-day mortality was considered as any death within the first month after 
103 the bloodstream infection onset. Persistent bacteremia was defined as the presence of positive 
104 blood cultures after 72 hours of appropriate antibiotic therapy. Septic embolism was 
105 determined as the presence of one or more diagnosed secondary foci as a result of bacterial 
106 spread through blood. Strains isolated from March 14th 2020 (Spain’s lockdown start) onwards 
107 were considered as pandemic clones. These data are available in Data S1. 
108 Strain isolation and manipulation. For each patient, only the first episode of S. aureus 
109 bacteremia was included in the analysis. Blood cultures were processed using the BACT/ALERT® 
110 automated system (bioMérieux) following standardized procedures. S. aureus strains were 
111 identified by using mass spectrometry technology (MALDI-TOF MS, Bruker). S. aureus single 
112 colonies were frozen and stored at -80°C. Antibiotic susceptibility testing was assessed by 
113 microdilution methodology according to standardized protocols using the MicroScan WalkAway 
114 — system" (Beckman Coulter). 
115 Ethics Statement: The Ethics Committee for Investigation with medicinal products (CElm) of the 
116 Parc Tauli University Hospital approved the implementation of this study (2023/5088, approved 
117 on 6 October 2023). The requirement for informed written consent was waived given the 
118 retrospective nature of the study. Patient identification was encoded, complying with the 


119 requirements of the Spanish Organic Law on Data Protection 15/1999. 
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120 


121 Whole-genome sequencing and de novo assembly 


122 DNA extraction. S. aureus isolates were grown on Columbia Agar with 5% Sheep Blood plates 
123 (bioMérieux) at 37°C. Total DNA was purified using the DNeasy Blood & Tissue Purification Kit 
124 (Qiagen) following the manufacturer’s instructions. Cell lysis was achieved by pre-incubating 
125 bacterial colonies in Phosphate-buffered saline (PBS) containing 100 ug/mL lysostaphin (Merck) 
126 at 372C for 30 minutes. DNA quality was assessed using a NanoDrop device (Thermo Fisher 
127 Scientific) and a Qubit® 2.0 fluorometer (Thermo Fisher Scientific). 

128 DNA sequencing. Libraries for sequencing were prepared with the Nextera XT DNA Sample 
129 Preparation Kit (Illumina). Whole-genome sequencing was performed using paired-end 
130 sequencing on an Illumina HiSeq 2500 and NovaSeq 600 sequencers available at the Genomics 
131 Unit of the Centre de Regulació Genòmica (CRG, Barcelona). Whole-genome sequencing data 
132 are available at National Center for Biotechnology Information (NCBI) under the accession 
133 number PRJNA1055690. 

134 Quality control. The quality of the raw sequencing reads was checked with FastQC (Andrews 
135 2010). Read pre-processing and filtering were performed with TrimGalore (Krueger 2012) by 
136 shaving the sequencing adapters, trimming the initial 20 bp poor-quality positions and using a 
137 Phred score >= 20 limit. 

138 De novo assembly. The trimmed paired-end reads were assembled de novo using shovill 
139 (Seemann 2023) and annotated with prokka (Seemann 2014) against the COG (Galperin et al. 
140 2021), HAMAP (Pedruzzi et al. 2015) and Pfam (Mistry et al. 2021) databases and using the S. 
141 aureus NCTC 8325 [NC_007795] proteome as the reference. CheckM2 was finally used for 


142 assessing the quality of assemblies (Chklovski et al. 2023). 


143 


144 Isolate typing and S. aureus annotation 


145 Isolate typing. Multilocus sequence typing (MLST) was carried out on the assembled scaffolds 
146 using the MLST software (Larsen et al. 2012). Sequence Types (ST) and Clonal Complexes (CC) 
147 were deduced from the S. aureus PubMLST database (Jolley et al. 2018) grouping the arcC, aroE, 
148 glpF, gmk, pta, tpi and ygiL gene types. spaTyper was used to assign spa types according to the 
149 Ridom Spa Server database guidelines (Sanchez-Herrero 2020), AgrVATE was used for agr typing 
150 (Raghuram et al. 2022) and SCCmec-types were obtained using SCCmecFinder (Kaya et al. 2018). 


151 Antibiotic resistance genes and virulence factors. Putative antibiotic resistance genes (ARG) 


152 and virulence factors (VF) were predicted on the assembled scaffolds with ABRicate (Seemann 
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153 2017) using the in-built NCBI AMRFinderPlus database (Feldgarden et al. 2019) and Virulence 
154 — Factor Database (VFDB) (Liu et al. 2019), respectively. Strains were considered as Panton- 
155 Valentine leukocidin (PVL) positive when both /ukF and lukS genes were predicted on the same 


156 genome. 


157 Plasmid prediction. Plasmid scaffolds were identified using MOB-recon and mobility was 
158 predicted based on the MOB-typer module (Robertson and Nash 2018). PlasmidFinder software 
159 — from Center for Genomic Epidemiology (CGE) was used for rep typing (Carattoli et al. 2014) and 


160 roary for identifying plasmid core genes (Page et al. 2015). 


161 Prophage prediction. Prophage prediction was performed as described recently for other S. 
162 aureus genomes (Sweet et al. 2023). Prophage regions were detected in assembled scaffolds 
163 using PhiSpy (Akhter et al. 2012) and CheckV was then used to remove host contamination and 
164 — to delimit phage boundaries (Nayfach et al. 2021). Prophage clustering was carried out with 
165 usearch (Edgar 2010) so that sequences showing 90% similarity along their 90% length were 
166 counted as the same. Integrase detection and typing was assessed with BLASTP (Altschul et al. 
167 1990) using previously reported integrase sequences (Goerke et al. 2009) as queries and limiting 
168 the e-value to le 0 and query coverage to >75%. Each putative integrase identified in the 
169 prophage regions was classified into an integrase group according to the query with the lowest 


170 e-value in the BLASTP analysis. 


171 Statistical methods 

172 Mann-Whitney U (MWU) test was used for continuous variables and Fisher’s exact test (FT) for 
173 categorical variables. For clonal associations, MWU was first applied to determine biased 
174 distributions of lineages when analyzing a specific genetic feature. If the lineage distributions 
175 differed significantly when compared to global frequencies, FT was independently applied to 
176 determine the association of each feature with specific lineages. Clonal associations were only 
177 computed for strains with known lineage. Always a p-value threshold of 0.05 was considered as 
178 significant. Statistical analyses were performed using custom Python scripts. 

179 

180 

181 RESULTS 


182 Cohort description. We identified 339 consecutive adult patients diagnosed with S. aureus 
183 bloodstream infections between 2014 and 2022 (Data S1). Most of them were male (69.3%) and 
184 — the mean age was 63.7 + 19.9 years. The calculated Modified Charlson Score was greater or 


185 equal than 4 in 64.0%. The most frequent clinical conditions detected were diabetes (31.6%), 
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186 the presence of a foreign device (25.4%), ischemic heart disease (17.7%), congestive heart 
187 failure (16.5%), solid neoplasm (15.3%), peripheral arteriopathy (12.7%) and chronic obstructive 
188 pulmonary disease (10.0%). Acquisition was nosocomial in 56.6%, while 27.4% were community- 
189 acquired and 16.0% healthcare related. Regarding the primary source of infection, vascular 
190 catheters (31.8%) and skin infections (18.0%) were predominant, followed by pneumonia (9.4%). 
191 As far as 30.7% episodes were of unknown origin (n=104). Antibiotic susceptibility tests indicated 


192 — that 13.9% infections were caused by methicillin-resistant S. aureus (MRSA) strains (Data S1). 


193 Regarding clinical outcomes, 93 cases were fatal (27.4%) and the median time to patient’s death 
194 was 20 (interquartile 25 (IQ25)-8, interquartile 75 (IQ75)-37) days (Data S1). Fatal bacteremia 
195 episodes showed to be significantly associated with age, modified Charlson, sepsis, congestive 
196 heart failure, immunosuppressive chemotherapy treatment, endocarditis and the presence of a 
197 cardiac device (Table S1). Of note, catheter-related bacteremia episodes were less lethal (p < 
198 0.05) (Table S1). On the other hand, persistent bacteremia was observed in 19.5% of episodes, 
199 with the median duration being 43 (1Q25=24, 1Q75=67.5) days (Data S1). Persistent infections 
200 were significantly associated with sepsis (Table $1). Pneumonia-originated bacteremia never 
201 resulted in persistent infections or septic emboli infections (p < 0.05) (Table S1). Septic emboli 
202 were diagnosed in 24.2% (Data S1) which were more prevalent in community-acquired 
203 infections (p < 0.05) (Table S1). Septic emboli infections were significantly associated with 


204 infective endocarditis episodes (Table S1). 


205 Isolate typing. The MLST software classified our S. aureus isolates in 24 different lineages (Figure 
206 1A, Data S1). The major identified lineages were CC30 (19.5%), CC5 (17.1%), CC45 (12.1%), ST398 
207 (8.8%), CC8 (8.8%), CC22 (8.0%) and CC15 (5.9%). 18 isolates could not be classified into a known 
208 sequence type. The distribution of lineages remained relatively stable over the years of study 


209 (Figure 1B). 


210 We found significant associations between specific lineages and clinical data. CC5 and CC8 
211 strains were more frequently recovered from patients with pneumonia-related bacteremia, 
212 while CC30 isolates were more prevalent among bacteremia from surgical site infections (p < 
213 0.05) (Figure S1, Table S2). Although we did not find any significant association, CC1 and CC8 
214 infections showed the highest mortality rates (45.4% and 40.0%, respectively). CC15, CC30, 
215 CC45, CC21 and ST398 lineages generated persistent infection with frequencies ranging from 
216 20.0 to 25.0%, while CC1 (9.1%) and CC22 (7.4%) were the lineages generating fewer persistent 


217 infections. The strains most frequently associated with septic emboli were CC15 and CC45 
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218 — (30.096 and 34.1%) while those that formed the least were CC1 and CC8 (9.1% and 10.0%) (Data 
219 S1). 


220 MRSA strains belonged more frequently to CC5 (42.5%) and CC8 (38.3%) lineages (p < 0.05) 
221 (Figure S1, Table S2) but methicillin resistance was also found in CC22 (8.5%), CC45 (4.2%) strains 
222 and unknown lineages (6.5%) (Data S1). The vast majority (95.7%) of SCCmec cassettes were 
223 grouped into type IV, except for two SCCmec cassettes classified as type V (Data S1). SCCmec 
224 type IVc was associated with CC5 and CC8 while SCCmec type IV was linked to CC5 and CC22 
225 clones (p < 0.05) (Figure S1, Table S2). In our collection, we did not find any association of MRSA 
226 strains with mortality (Matthews Correlation Coefficient (MCC) = 0.06). 


227 Our analysis predicted that 13 isolates from our collection were agr-negative. Among the agr- 
228 positive strains, 48.5% were classified as type |, 26.1% type II, 22.1% type Ill and 3.3% type IV 
229 (Data S1). We observed a biased distribution of lineages regarding the agr types II and III (p < 
230 — 0.05) (Figure S1). Specifically, we found that type ll agr was associated with CC5 and CC15 
231 isolates, while type Ill was almost exclusive of CC30 clones (p < 0.05) (Table S2). Finally, our 
232 typing analysis revealed 123 different spa types (Data S1), being t012 significantly associated 
233 with CC30 isolates, t002 and t067 with CC5, t008 with CC8, t084 with CC15 and t1451 with ST398 
234 (Figure S1, Table S2). 


2014 2015 2016 2017 2018 2019 2020 2021 2022 


m CC30 MB CC5 E CC45 mST398 mCC8 mCC22 mCC15 mCC1 MmCC121 m Other 


235 


236 Figure 1. (A) Pie chart depicting global S. aureus lineage frequencies in our collection. (B) 
237 Graphical representation of S. aureus lineage frequencies over years. Lineages with less than 5 


238 — representatives are collected under the umbrella term "Other". 


239 Antibiotic resistance genes. We identified homologs to 34 different ARG families conferring 


240 resistance to clinically relevant antibiotics (Figure 2, Data S1). Homologs of the mecA gene were 
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241 found in 47 strains (13.9%) and they were associated primarily with CC5 and CC8 strains (p < 
242 — 0.05) (Figure S1, Table S2). 


243 Homologs of genes related to tetracycline, B-lactams and fosfomycin resistance were the most 
244 prevalent. Homologs of tet38 were identified in all isolates, of complete blal-blaR-blaZ clusters 
245 in 279 strains (82.3%) and of fosB in 206 (60.2%). We found that the blaZpc: variant was mainly 
246 associated with strains belonging to the CC30 lineage (p < 0.05) (Figure S1, Table S2). All isolates 
247 belonging to CC5, CC8, CC15, CC30 and CC121 lineages encoded a fosB gene homolog. In 
248 contrast, CC22, CC45 and ST398 clones did not present any fosfomycin resistance gene 
249 homologs (Table S2). 


250 In addition, we identified homologs to three different groups of genes conferring macrolide 
251 resistance: several erm variants encoding 23S rRNA methylases (19.5%), the ABC-F type 
252 ribosomal protection encoding gene msrA (10.0%) and the phosphotransferase gene mphC 
253 (9.4%). We observed a strong correlation between homologs of mphC and msrA (MCC = 0.97), 
254 which were significatively associated with CC5 and CC8 lineages (Figure S1, Table S2). In addition, 
255 almost all strains carrying ermT belonged to the ST398 sequence type (p < 0.05) (Figure S1, Table 
256 S2). Homologs of ermA were associated with CC1 and CC5 while ermC homologs were linked to 
257 —CC8 clones (p < 0.05) (Figure S1, Table S2). We also detected a biased distribution of the 
258 aminoglycoside nucleotidyltransferase genes ant(4’)-la (9.1%) and ant(9)-la (5.9%) as well as of 
259 — the aminoglycoside phosphotransferase genes aph(3')-Illa (7.4%) and aph(2’’)-Ih (4.4%) (Figure 
260 = $1). Specifically, ant(4’)-la gene homologs were significatively associated with CC5 and CC8 
261 strains, ant(9)-la with CC1 and CC5, aph(3'J-Illa with CC5 and CC8, and aph(2")-Ih with ST398 
262 (Table S2). Homologs of the bleomycin resistance gene bleO (7.4%) showed a positive correlation 
263 with ant(4’)-la (MCC = 0.89) and were also linked to CC5 and CC8 clones (p < 0.05) (Figure S1, 
264 Table S2). Finally, homologs of the mupA gene (3.8%), conferring mupirocin resistance, were 


265 —significatively associated with CC30 isolates (Figure S1, Table S2). 
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267 Figure 2. Clustered heatmap representing the frequency of abricate-predicted ARGs for all S. 
268 aureus lineages. Cells are colored from pale orange (0%) to black (100%). Samples were clustered 
269 and visualized using the Seaborn clustermap function. Isolates with unknown lineage were not 


270 included. 


271 Virulence factors. The presence of VF contributing to S. aureus pathogenicity and invasiveness 
272 was predicted (Figure 3, Data S1). All isolates from our collection carried homologs of the cap 
273 — gene cluster responsible for the synthesis of capsular polysaccharide. Capsular polysaccharide 
274 serotype 8 was significatively more prevalent in CC15, CC30, CC45 and CC121 lineages (Figure 
275 S1, Table S3) while serotype 5 were mainly found at CC5, CC8, CC22 and ST398 clones. Homologs 
276 ofthe Immune Evasion Cluster (IEC) (sak, chp and scn) and sbi, also involved in immune system 
277 evasion, were present in 73.1 - 99.4% of the isolates. The sak gene was not found in ST398 
278 isolates. Similarly, 99.7% of the isolates carried the isdA-G gene homologs, which confer 
279 resistance to killing by human lactoferrin. The intracellular adhesion (ica) locus, involved in cell- 
280 cell adherence and biofilm formation, was detected in all strains. Homologs of genes coding for 
281 several Microbial Surface Components Recognizing Adhesive Matrix Molecules (MSCRAMMs) 
282 were also detected, although to different extents: ebp (95.3%), eap/map (93.5%), sdrC (92.0%), 
283 —sdrE (87.6%), clfA (86.1%), fnbA (74.6%), cIfB (68.4%), fnbB (62.2%), sdrD (53.7%), vWbp (47.8%) 
284 and cna (15.9%). Specifically, we found that the secreted von Willebrand factor-binding protein 
285 (vWbp) was more prevalent in CC1, CC8, CC30, CC45 and CC121 lineages while the cna gene was 
286 associated with CC30 clones (p < 0.05) (Figure S1, Table S2). 


10 
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287 = Almost all isolates carried homologs of the genes encoding the a-, B- and ô- hemolysins (hla, hlb, 
288 —hIgABC and hid). Genes responsible for the formation of type VII secretion systems (T7SS), esa, 
289 ess and esx, were also widespread in our collection. Enterotoxin-encoding genes (sea, seb, sec, 
290 sed, seh, selk, sell and selq) showed markedly different frequencies ranging from 31.0% (sea) to 
291 2.9% (selk and selq) and significatively biased distributions among different lineages (Figure S1). 
292 Specifically, CC1 strains were associated with seb, seh, selk and selq; CC5 with sea and sed; CC8 
293 with sea, sell and selq; CC30 with sea and seh; CC45 with sea, sec and sell, and CC121 with seb 
294 (p< 0.05) (Table S2). Of note, sec and sell genes were totally correlated (MCC = 1.0). We also 
295 detected the presence of the toxic shock syndrome toxin-1 (tsst-1) gene homologs in 56 isolates 
296 — (16.596), which were linked to CC22 and CC30 (p < 0.05) (Figure S1, Table S2). In addition, 11 
297 isolates were predicted to encode the PVL genes (/ukF and /ukS) with a linkage to CC8 strains 
298 (Figure S1, Table S2). Finally, we also found homologs of the toxin-encoding genes eta (2.1%) 
299 — and etb (0.9%) in our strains. 
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302 Figure 3. Clustered heatmap representing the frequency of abricate-predicted VFs for all S. 
303 aureus lineages. Cells are colored from pale orange (0%) to black (100%). Samples were clustered 
304 and visualized using the Seaborn clustermap function. Isolates with unknown lineage were not 


305 included. 


306 Plasmids. We identified 51 different plasmids distributed in 256 isolates (75.5%) (Figure 4, Data 
307 S1). Among them, the vast majority (52.8%) presented a single plasmid, but we also found strains 


308 with two (18.3%), three (3.2%) or four (1.2%) plasmids. 
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309 = The core proteome of the most prevalent plasmid, AA410 (n=72; rep5a/rep16/rep19), showed 
310 heavy metal resistance (arsB, arsC, cadA, cadC, copB and mco) and lactococcin synthesis genes 
311 (Table S3). This plasmid was linked to CC30 strains (p < 0.05) (Figure $1, Table S2). Plasmid AA411 
312 (n=56; rep5a/rep16), harbored the blal-blaR-blaZ cluster for penicillin resistance, a cadmium 
313 resistance cluster (cadC, cadD) and lactococcin-synthesis genes (Table S3). In this case, AA411 
314 plasmids were significatively associated with CC1, CC15 and CC45 (Figure S1, Table S2). Plasmid 
315 AA840 (n=33; rep16/rep19/rep21), encoding cadA gene for cadmium resistance and several 
316 ARGs (blal-blaR-blaZ, aph(3')-Illa, mphC, msrA and sat4) (Table S3), was more prevalent among 
317 —CC5 and CC8 isolates (p < 0.05) (Figure S1, Table S2) and showed a positive correlation with 
318 MRSA strains (MCC = 0.73). Also significatively associated with these lineages (Figure S1, Table 
319 S2) was the AA851 plasmid (n=22, rep22) which harbored cadmium resistance genes (cadC, 
320 cadD), the blal-blaR-blaZ cluster and genes for lactococcin synthesis (Table S3). On the other 
321 hand, the plasmid AB641 (n=25; rep13), which was almost exclusively detected in isolates from 
322 the ST398 lineage (p < 0.05) (Figure S1, Table S2), carried two cadmium resistance genes (cadC, 
323 cadD) and the macrolide-resistance gene ermT (Table S3). Finally, plasmid AC670 (n=19; 
324 —rep7c/rep16/rep19/rep22) encoded genes conferring antimicrobial resistance (ant(4')-la, bleO); 
325 and plasmid AA848 (n=15; rep20) carried several ARGs (blal-blaR-blaZ, mphC, msrA and tetk), 
326 two cadmium resistance loci (cadC, cadD) and the staphylococcal enterotoxin D-encoding (sed) 
327 gene (Table S3). AC670 plasmids were associated with CC5 and CC8 clones while AA848 plasmids 
328 were linked to CC5 clones (p < 0.05) (Figure S1, Table S2). 


329 We also predicted other less-abundant plasmids, carrying other genes conferring resistance to 
330 aminoglycosides (ant(6)-la, aph(2")-Ih), macrolides (ermA, ermB, ermC, msrA, mphA), 
331 tetracycline (tetR), chloramphenicol (catA7), lincosamides (InuA) and trimethoprim (dfrG). The 


332 exfoliative toxin B encoding gene (etb) was also annotated in three plasmids (Table S3). 


333 Overall, our results showed that ant(6)-la, aph(3')-llla, catA7, ermC, InuA, mphC, msrA and sat4 
334 genes were exclusively plasmid-borne in our collection (Figure S2). Likewise, mupA (92.3%), 
335 ant(4')-la (90.3%), bleO (88.0%), ermT (80.6%) and tetR (77.8%) genes were also much more 
336 prevalent in plasmids than in chromosomes. Similarly, around 50% of the blal-blaR-blaZ cluster 
337 was detected in plasmids. In contrast, several abundant genes such as ant(9)-la, ermA, fosB, 
338 mecA or tet38 were always chromosomally encoded (Figure S2). Finally, the virulence genes sed 


339 and etb were found to be exclusively plasmid-borne. 
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340 


341 Figure 4. Schematic representation of the most abundant plasmids in our collection. Gene 
342 coloring depicts ARGs (gray), VF (orange), heavy metal resistance genes (purple) and lactococcin 
343 synthesis genes (white). Rep types and lineage associations are also included. Only plasmids 


344 harbored in more than 15 isolates were included. 


345 Prophages. We predicted 369 different prophage clusters in 336 isolates (97.6%) (Data S1). The 
346 majority presented one (28.6%) or two (28.9%) prophages, but we also found strains with three 
347 — (21.5%), four (12.1%), five (7.4%), six (0.9%) or even seven (0.6%) prophage sequences. 


348 Prophage clusters from the integrase group 3 (Sa3) were predominant in our collection (n=43). 
349 Some of them (n=15) carried homologs of complete IEC (scn, chp, sak). Enterotoxin A encoding 
350 gene (sea) was predicted in six of them and tsst-1 was found in one cluster (Table S4). Of note, 
351 different Sa3-prophage clusters (1, 15, 29 and 55) showed a biased distribution in our collection 
352 (Figure S1). Specifically, cluster 1 was associated with CC30, cluster 15 with ST398, cluster 29 
353 — with CC5 and, finally, cluster 55 with CC30 and ST398 strains (p < 0.05) (Table S2). 


354 We further identified some prophage clusters (n=26) of the integrase group 2 (Sa2), two of them 
355 encoding sea and tsst-1 genes. We also predicted 24 prophage clusters from integrase group 1 
356 (Sai); 22 from integrase group 5 (Sa5), with two harboring sec and sell genes; and 17 from 
357 integrase group 6 (Sa6). Prophages from integrase group 7 (Sa7) and 9 (Sa9) were also identified 
358 but to a lesser extent. Finally, 221 prophage clusters did not present any integrase gene. Among 


359 those, 14 clusters carried a complete IEC and 10 harbored the sea gene (Table S4). 


360 No ARG were detected in the predicted prophage genomes. 
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361 Pandemic-associated features. We identified 139 and 200 consecutive adult patients diagnosed 
362 with bloodstream infections by S. aureus before and after March 14th 2020, respectively (Data 
363 S1). In the pre-pandemic years, we detected 27.8 episodes of S. aureus bacteremia per year in 
364 — our hospital rising to 66.7 during the pandemic. In such years, we observed an increase of 
365 nosocomial infections, from 69 episodes in the pre-pandemic years to 119 during the COVID-19 
366 pandemic. Specifically, cases of S. aureus peripheral catheter-related bacteremia were 
367 significantly more prevalent and we also found 2.5 times more episodes of bacteremia 
368 originating from surgical sites (Table S5). We also identified a significant increase of MRSA 


369 infections (p < 0.05) and we diagnosed 2.8 times more endocarditis cases (Table S5). 


370 Statistical analysis on pre- and pandemic strains revealed significant genetic differences 
371 between these two S. aureus populations. We observed a strong increase of strains carrying 
372 genes encoding the clumping factors A (c/ifA, p = 2.1e S) and B (c/fB, p = 2.3e°°) (Table S5). 
373 Clumping factors A and B are MSCRAMMs that bind fibrinogen, which is found at high levels 
374 upon SARS-CoV-2 infection (O’Brien et al. 2002; Kangro et al. 2022). In contrast, we found a 
375 significant reduction of strains encoding the sdrE gene, which is also classified as MSCRAMM 
376 (Table S5). On the other hand, we noticed a rise in several antibiotic resistance genes during 
377 those years. Specifically, the macrolide-resistance msrA and mphC genes as well as the 
378 methicillin-resistance mecA gene showed a significant increase during the pandemic years 
379 (Table S5). Of note, macrolides were extensively used as accessory therapy for COVID-19 in 


380 patients (Sterenczak et al. 2020; Grau et al. 2021). 
381 
382 DISCUSSION 


383 Global insights into S. aureus clinics, antibiotic resistance and virulence 

384 Genomics enables unraveling the genetic makeup of pathogens, providing insights into virulence 
385 factors, antibiotic resistance mechanisms and the overall pathogenic potential of bacteria. This 
386 information is crucial for the development of targeted diagnostics, treatment strategies and 
387 preventive measures. The use of genomics also allows tracking the spread of specific strains and 
388 identifying emerging threats, contributing to the surveillance and control of infections. Here, 
389 using a collection of 339 S. aureus isolates from patients with bacteremia, we unveiled important 
390 insights into the clinical and genomic features of these bacterial pathogens (Table 1). 

391 In our hospital, we identified S. aureus lineages mirroring those previously reported in Spain 
392 —(Pérez-Montarelo et al. 2018) and in other countries (Fowler et al. 2007; Park et al. 2017; 


393 Campbell et al. 2022) (Figure 1), suggesting a shared genetic heritage and potentially common 
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394 evolutionary pathways. Pérez-Montarelo and colleagues analyzed strains isolated across various 
395 Spanish hospitals spanning 2002 to 2017 (Pérez-Montarelo et al. 2018). Our study, conducted 
396 between 2014 and 2022, complements this report and broadens its scope into the pandemic 
397 period, offering compelling observations on the evolutionary trends of S. aureus clones. 
398 Specifically, we observed an increase of ST398 lineages, heightening from 2.4% (Pérez- 
399 Montarelo et al. 2018) to a mean prevalence of 8.8% in our isolates. Importantly, the prevalence 
400 ofthis lineage in our collection ranged from 0% (2014 and 2015) to 15% (2022). These results 
401 highlight the spread of ST398 lineage in hospital settings, confirming the need to focus on the 
402 mechanisms involved in the epidemiological success of ST398 in humans, as previously 
403 suggested (Uhlemann et al. 2017; Mama et al. 2021; van der Mee-Marquet et al. 2022) (Figure 
4041). 

405 Despite possible associations between genotypes and certain complications, the role of S. 
406 daureus genetic factors in the evolution of bacteremia is still poorly understood. CC1, CC5, CC8, 
407 CC22 and CC30 lineages have been associated with increased risk of complications (Fowler et al. 
408 2007; Miller et al. 2012). In contrast, more recent studies did not find any apparent association 
409 between clonality and poorer outcomes (Chong et al. 2013; Gasch et al. 2014; Fernandez- 
410 Hidalgo et al. 2018). Our study is in line with the latter since we did not find significant 
411 associations between the different isolates from our collection and clinical outcome. However, 
412 there is a tendency for the CC1 and CC8 strains to cause fatal infections. On the other hand, we 
413 identified associations between the different lineages and the focus of infection. CCS and CC8 
414 were correlated with pneumonia-originated bacteriemia, as deduced from previous reports 
415 (Antonelli et al. 2019), and CC30 with surgical site infections (Table 1, Table S2). 

416 MRSA poses a significant public health concern due to its ability to cause infections that are 
417 challenging to treat, showing higher mortality than MSSA infections (Munckhof et al. 2008). In 
418 Spain, the prevalence of methicillin resistance among S. aureus isolates is around 20-25% (Pérez- 
419 Montarelo et al. 2018; Vazquez-Sanchez et al. 2022), while the incidence of MRSA infections in 
420 our hospital showed to be 13.9% (Figure 2). Of note, we did not find associations with mortality 
421 and MRSA strains. In accordance with some reports, most of our MRSA isolates belonged to CC5 
422 or CC8 strains (Dupper et al. 2019). These lineages were also the ones accumulating more 
423 antibiotic resistance genes (Figure 2, Table 1), as previously stated (Smith et al. 2021). Genes 
424 predicted to confer resistance to tetracycline and B-lactams were the most prevalent in our 
425 isolates (Figure 2). However, the results for tetracycline did not align well with the patterns 


426 observed in vitro by broth microdilution tests. The tet38 gene, which was present in almost all 
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427 strains, is known to confer tetracycline resistance when overexpressed (Chen and Hooper 2018), 


428 but its mere presence does not necessarily indicate resistance to this class of antibiotics. 


429 S. aureus is known for its pathogenic potential, primarily attributed to a diverse set of virulence 
430 factors. As previously elucidated, the agr groups and capsular polysaccharide serotypes were 
431 associated with S. aureus clonality (Figure 3, Table 1) (Monecke et al. 2011). Several genes 
432 involved in immune system evasion (adsA, aur, sbi, ssl, sspA, sspB), lactoferrin resistance (isdA- 
433 G), adhesion (ica), T7SS formation (esa, ess, esx), some MSCRAMMs (ebp, eap/map, sdrC) and 
434 hemolysins were found in virtually all strains (Figure 3), similar to previous reports (Rasmussen 
435 et al. 2013). IECis an immune response cluster harboring chp, scn and sak genes which is typically 
436 encoded by B-hemolysin converting bacteriophages (Sieber et al. 2020). This is partially 
437 paralleled in our isolates, in which we predicted 61.7% of complete IEC clusters in prophage 
438 regions (Table S4). In contrast to many other bacterial pathogens, S. aureus produces a wide 
439 array of toxins (Cheung et al. 2021). In our hospital, besides hemolysins, the most abundant 
440 predicted toxins were the exotoxins SEA (31.0%) and TSST-1 (16.5%), the latter being the major 
441 cause of toxic shock syndrome (Figure 3). SEA is the most common staphylococcal enterotoxin 
442 (Pinchuk et al. 2010) and is also reported to be encoded by a prophage (Betley and Mekalanos 
443 1985). In our isolates, only 53.8% of sea genes were predicted in prophage regions (Table S4), 
444 hereby indicating that we could be underestimating prophage regions. Finally, we found 11 PVL- 
445 positive strains (3.2%) (Figure 3). Panton—Valentine leukocidin is a pore-forming leukotoxin 
446 produced by 2-3% of S. aureus isolates (Wojcik-Bojek et al. 2022) which has been traditionally 
447 used as an indicator of community-acquired MRSA strains (Fey et al. 2003). However, some 
448 recent studies found that PVL may not be a reliable marker because it is being reported in 
449 nosocomial MRSA infections (Nakaminami et al. 2020; Shohayeb et al. 2023). Paralleling these 
450 reports, we found 4 non-community MRSA harboring PVL suggesting a possible spread of these 


451 strains in hospitals. 
452 
453 
454 
455 
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Table 1. Summary of the significant associations between the most abundant S. aureus lineages 


with all the clinical and genetic data. 


Lineage Clinical Gene Antibiotic | Virulence | Plasmids Prophages 
features types Resistance — Factors 
Genes 
cc1 ant(9)-la, ermA seb, seh, sell, | AA411 
selq 
CC5 MRSA, pneumonia- | SCCmec IVc, agr | ant(4’)-la, ant(9)- | cap5, sea, sed AA840, AA851, — cluster_29 (Sa3) 
originated bacteremia | type ll, spa t002, la, aph(3’)-Illa, AC670, AA848 
spa t067 bleO, ermA, fosB, 
mecA, mphC, 
msrA 
ccs MRSA, pneumonia- | SCCmec IVc, spa | ant(4’)-la, cap5, sea, selk, | AA840, AA851, 
originated bacteremia | t008 aph(3')-llla, bleO, | selq, PVL AC670 
ermG, fosB, 
mecA, mphC, 
msrA 
CC15 agr type ll, spa | fosB cap8 AA411 
t084 
cc22 SCCmec IV cap5, tsst-1 
cc30 Bacteriemia from | agr type Ill, spa | blaZPC1, fosB, | cap8, sea, seh, | AA410 cluster_1 (Sa3), 
surgical site infections | t012, mupA tsst-1 cluster_55 (Sa3) 
CC45 cap8, sea, sec, | AA411 
sell 
CC121 fosB cap8, seb 
ST398 spat1451 aph(2”)-Ih, ermT cap5 AB641 cluster 15 (Sa3), 


cluster 55 (Sa3) 


Changes in molecular characteristics of pandemic S. aureus isolates 


Understanding the interplay betvveen viral and bacterial infections, including the prevalence and 
impact of S. aureus in COVID-19 cases, is important for effective patient care and public health 
measures. Several studies underscored a higher incidence of bacteremia by S. aureus at the 
COVID-19 pandemic (Falces-Romero et al. 2023). In keeping with this, the episodes of S. aureus 
bacteremia per year doubled in our hospital during the pandemic. Concomitantly, we found a 
rise in nosocomial infections, particularly those originated from catheters and surgical sites. The 
increase of catheter-related bacteremia during the COVID-19 epidemic has been previously 


documented by our group (Gasch et al. 2022). 


Unlike other countries, we did not detect emerging S. aureus clones during the pandemic (Gu et 
al. 2023). However, we noticed a significant increase in strains carrying the c/fA and c/fB genes, 
a reduction of strains encoding sdrE, along with a rise in macrolide and methicillin antibiotic 


resistance genes (Table S5). 
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474 The cIfA and cIfB genes are MSCRAMMs encoding the fibrinogen-binding Clumping Factors A 
475 and B, respectively. Fibrinogen is one of the most prevalent coagulation proteins in blood 
476 (Kangro et al. 2022). CIfA is the major virulence factor responsible for S. aureus clumping in blood 
477 plasma (O’Brien et al. 2002) and CIfB interacts with cytokeratin 10 and loricrin facilitating S. 
478 aureus skin infection (Lacey et al. 2019). Both CIfA and CIfB have been reported to be involved 
479 in bacterial and platelet aggregation (O’Brien et al. 2002). Interestingly, several reports 
480 documented an increase of fibrinogen production by SARS-CoV-2 virus (Bouck et al. 2021; 
481 Kangro et al. 2022). Therefore, our results suggest that the presence of c/fA and c/fB genes 
482 provided an advantage to S. aureus in the context of SARS-CoV-2 co-infection by enhancing 
483 fibrinogen binding. In consequence, S. aureus bloodstream isolates carrying the c/fA and c/fB 
484 genes were positively selected between 2020 and 2022. Conversely, we found a reduced 
485 prevalence of the sdrE gene, which codes for a MSCRAMM that prevents S. aureus phagocytosis 
486 by sequestering the human complement factor H (CFH) on the surface of bacterial cells (Zhang 
487 et al. 2017, Wojcik-Bojek et al. 2022). This is remarkable because it has been shown that the 
488 SARS-CoV-2 spike protein blocks CFH leading to a complement dysregulation on the human cell 
489 surface (Yu et al. 2022). Therefore, the functional redundancy of S. aureus SdrE and the SARS- 
490 CoV-2 spike protein could explain the negative selection of sdrE. In the same line, it is important 
491 — to note that S. aureus CflA has also been found to have an anti-phagocytic effect (Higgins et al. 


492 2006; Hair et al. 2010). 


493 Regarding S. aureus antibiotic resistance, some articles have emphasized the role of macrolides 
494 as adjunctive therapy for COVID-19 in certain patients (Sterenczak et al. 2020) and the use of 
495 such antibiotics could allow selecting resistant strains. Macrolide resistance in staphylococci is 
496 mainly associated with mobile ARGs which often correlates with methicillin resistance (Serra et 
497 al. 2023), as also observed in our collection (Figure 2). These findings emphasize the evolving 


498 nature and impact of S. aureus infections in the COVID-19 pandemic era. 
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523 FIGURES AND TABLES 


524 Figure 1. (A) Pie chart depicting global S. aureus lineage frequencies in our collection. (B) 
525 Graphical representation of S. aureus lineage frequencies over years. Lineages with less than 5 


526 representatives are collected under the umbrella term “Other”. 


527 Figure 2. Clustered heatmap representing the frequency of abricate-predicted ARGs for all S. 
528 aureus lineages. Cells are colored from pale orange (0%) to black (100%). Samples were clustered 
529 and visualized using the Seaborn clustermap function. Isolates with unknown lineage were not 


530 included. 


531 Figure 3. Clustered heatmap representing the frequency of abricate-predicted VFs for all S. 
532 aureus lineages. Cells are colored from pale orange (0%) to black (100%). Samples were clustered 
533 and visualized using the Seaborn clustermap function. Isolates with unknown lineage were not 


534 included. 
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535 Figure 4. Schematic representation of the most abundant plasmids in our collection. Gene 
536 coloring depicts ARGs (gray), VF (orange), heavy metal resistance genes (purple) and lactococcin 
537 synthesis genes (white). Rep types and lineage associations are also included. Only plasmids 


538 harbored in more than 15 isolates were included. 


539 Table 1. Summary of the significant associations between the most abundant S. aureus lineages 


540 with all the clinical and genetic data. 
541 
542 SUPPLEMENTARY MATERIAL 


543 Data S1. JSON-formatted file including the clinical metadata and predicted features for all S. 


544 aureus Strains included in this study. 


545 Figure S1. Pie plots depicting lineage frequencies for all clinical and genetic features. The MWU 


546 p-values obtained when comparing to global lineage frequency are included. 


547 Figure S2. Heatmap representing the frequency of chromosome- or plasmid-encoded genes for 
548 all ARGs in S. aureus lineages. Cells are colored from blue (all instances are plasmid-borne) to 


549 red (all instances are in the chromosome), with white indicating absence of ARG. 


550 Table S1. Results for statistical analysis of clinical data. In each case, the computed p-value and 


551 selected statistical method are included. 


552 Table S2. Results of Fisher’s exact test (FT) for analyzing associations between specific S. aureus 
553 lineages and clinical or genetic data. FT was only computed for features with biased distributions 


554 — in MWU analysis (Figure S1). 


555 Table S3. List of predicted plasmids including number of instances, the roary-core proteome, 


556 and the abricate-predicted antibiotic resistance genes and virulence factors. 


557 Table S4. List of predicted prophage clusters including number of instances, integrase type, 
558 presence of IEC and the abricate-predicted virulence factors. All phage instances of the same 


559 cluster not necessarily harbored all the features. 


560 Table S5. Results of Fisher’s exact test for comparing the clinical and genetic features of pre- and 


561 pandemic bacteremia episodes. 
562 
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